/* 
 * File:   Gibbs.h
 * Author: hadrien
 *
 * Created on 15 avril 2012, 23:57
 */

#ifndef GIBBS_H
#define	GIBBS_H
#include <cstdlib>
#include "GaussianProcess.h"
#include <math.h>
#include <stdio.h>
#include <iostream>

template<typename ENERGYCALCULATOR, typename GIBBSPARAMETERS>
class Gibbs {
public:

    typedef ENERGYCALCULATOR energy_t;
    typedef typename energy_t::Params Params;
    typedef typename Params::features_vector_t features_vector_t;
    typedef typename Params::label_t label_t;
    typedef GIBBSPARAMETERS GibbsParams;
    typedef GaussianProcess<Params> model_t;

    void operator()(const model_t& model, features_vector_t const * const * observations, label_t** labels) const {
        int i, j;
        double *Ek; // array to store local energies
        label_t s;
        double summa_deltaE, E;
        double sumE;
        double z;
        double r;

        energy_t energy(&model, observations);

        //Initialisation to ML
        for (int i = 0; i < Params::getWidth(); ++i) {
            for (int j = 0; j < Params::getHeight(); ++j) {
                labels[i][j] = energy.getMaximumLikelyhood(i, j);
            }
        }

        Ek = new double[Params::getL()];

        int T = GibbsParams::getInitialTemperature();
        double E_old = energy.globalEnergy(labels);
        int it = 0;

        do {
            summa_deltaE = 0.0;
            for (i = 0; i < Params::getWidth(); ++i)
                for (j = 0; j < Params::getHeight(); ++j) {
                    sumE = 0.0;
                    for (s = 0; s < Params::getL(); ++s) {
                        Ek[s] = exp(-energy.localEnergy(i, j, s, labels) / T);
                        sumE += Ek[s];
                    }
                    r = ((double) rand()) / (double) RAND_MAX; // r is a uniform random number
                    z = 0.0;
                    for (s = 0; s < Params::getL(); ++s) {
                        z += Ek[s] / sumE;
                        if (z > r) // choose new label with probabilty exp(-U/T).
                        {
                            labels[i][j] = s;
                            break;
                        }
                    }
                }
            E = energy.globalEnergy(labels);
            summa_deltaE += fabs(E_old - E);
            E_old = E;
            T *= GibbsParams::getDecreaseTemperature(); // decrease temperature
            it++;
            printf("\nEnergy: %f\n", E);
        } while (summa_deltaE > GibbsParams::getStopCriterion()); // stop when energy change is small
        printf("\nNumber of iteration : %d, Energy: %f\n", it, E);
        delete[] Ek;
    }
private:

};

#endif	/* GIBBS_H */

